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Room temperature ionic liquids (RTILs) are solvent with unusual properties, which are difficult 
to characterize experimentally because of their intrinsic complexity (large number of atoms, strong 
Coulomb interactions). Molecular simulations have therefore been essential in our understanding of 
these systems. Depending on the target property and on the necessity to account for fine details of 
the molecular structure of the ions, a large range of simulation techniques are available. Here I focus 
on classical molecular dynamics, in which the level of complexity of the simulation, and therefore 
the computational cost, mostly depends on the force field. Depending on the representation of the 
ions, these are either classified as all-atom or coarse-grained. In addition, all-atom force fields may 
account for polarization effects if necessary. The most widely used methods for RTILs are described 
together with their main achievements and limitations. 


I. INTRODUCTION 


RTILs are a particular class of solvents in which all 
the species are ionic. They have started to be widely 
investigated much more recently than conventional sol¬ 
vents such as water or organic liquids. As a consequence, 
unlike the latter which have generally been studied by 
theoretical approaches long after many experiments had 
been performed, RTILs are an interesting testing field 
for the predictive power of molecular simulations. In ad¬ 
dition, most of the liquid state theories do not hold in 
these media due to the predominance of Coulombic in¬ 
teractions. A deep undestanding of the RTILs properties 
is nevertheless necessary for their use in a variety of ap¬ 
plications. Indeed, although their industrial use is cur¬ 
rently limited to the role of solvent for organic reactions, 
their good stability and their wide electrochemical win¬ 
dow makes them excellent electrolytes, e.g. for energy 
storage devicesThey are also able to stabilize new 
products such as nanoparticles with differentproperties 
than the ones produced in aqueous solvents.There are 
thousands of possible cat ions/anions combinations, and 
the knowledge of the structure and transport properties, 
in the bulk or at interfaces, is mandatory for choosing the 
most appropriate one. Therefore, many studies involve 
the use of molecular simulations in order to understand 
and predict such properties. 

Although quantum chemistry methods are very use¬ 
ful for characterizing the nature of the interactions in 

RTILsjStiSl 

these techniques are limited to small gas 
phase systems and they are not yet appropriate for the 
sampling of liquid properties. Most of the literature 
therefore involves molecular dynamics (MD) simulations. 
Interestingly, the outburst of ab initio molecular dynam¬ 
ics (AIMD), in which the atomic forces are calculated 
at each time step of the simulation using an electronic 
structure Density Functional Theory (DFT) calculation, 
has not been very visible in RTILs. There are two main 


explanations for this: Firstly, these liquids are strongly 
short-ranged ordered, each ion being surrounded by suc¬ 
cessive shells of alternating charged species. A relatively 
large number of ion pairs therefore have to be included 
in the simulation cells in order to avoid finite-size effects. 
Secondly, each molecule contains a large number of atoms 
and thus of electrons. The computational cost is there¬ 
fore particularly large, compared to volvent-based elec¬ 
trolytes. The few AIMD studies have therefore focussed 
on the calculation of properties which are out of reach of 
classical MD, such as chemical react ivit}Wl^ or vibra¬ 
tional propertieJ^. 
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FIG. 1: Non-exhaustive summary of the system sizes (number 
of ion pairs in the simulation cell) and the trajectory length 
simulated in a series of recent MD works. “AI-MD” stands 
for Ab initio MD, while “Class-MD” and “Pol-MD” respec¬ 
tively stand for classical MD involving non polarizable and 
polarizable force fields. Data extracted from reference^^^^^. 

The separation of time scales and length scales which 
can be studied by AIMD and classical MD is illustrated 




on figure which reports data from a series of articles 
(from the period 2012-2014) for the number of ion pairs 
in the simulation cell with respect to the length of the 
trajectories in the production runs. It is not an exhaus¬ 
tive set of data and there are many biases which have to 
be kept in mind: To cite only one, AIMD is more sen¬ 
sitive to the number of electrons than to the number of 
ions. Nevertheless, it clearly shows that there are two or¬ 
ders of magnitude of difference in the size and time that 
can be studied with the two methods. Despite the de¬ 
velopment of very efficient codes, it is very unlikely that 
AIMD will be able to bridge this gap before long, so that 
classical MD will remain in future years the method of 
choice for many problems involving RTILs and requiring 
long simulation times. 

Although they share a similar methodology, i.e. the 
iterative integration of Newton’s equation of motion over 
several millions of femtosecond timesteps, there exist 
many flavors of classical MD simulations. They are gen¬ 
erally classified according to the nature of the interacting 
bodies and on the analytical expression for the interac¬ 
tion potential (or force held) which acts between them. 
For the former, the most general way to handle molecu¬ 
lar species is to assign one interaction site to each atom 
(all atom force fields). This approach may be simpli¬ 
fied by coarse-graining, for example by gathering the hy¬ 
drogen atoms with the atoms to which they are linked 
(united atoms force fields), or even further by including 
a larger number of atoms into one grain (coarse-grained 
force fields). There is an even larger range of choices for 
the analytical expression of the interaction potential. In 
this article, the most widely used methods for RTILs will 
be described. The following key points will be addressed: 
In each case, the general shape of the force held will be 
provided, together with some general methods for obtain¬ 
ing the corresponding parameters. Then the capabilities 
and limitations of each method will be discussed. 


II. NON-POLARIZABLE ALL ATOM FORCE 
FIELDS 


take the following forms: 
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The sums respectively run over the full set of instanta¬ 
neous bonds angles Oijk and dihedrals <pijki which 
must be defined prior to the simulations (i.e. the bonds 
are not allowed to break or form during the simulation). 
The quantities 0^-^ and are parameters 

of the force held. 

The non-bonded term is in principle composed of three 
contributions 


-S'non—bonded — A^repulsion T A^dispersion T -^^electrostatic (5) 

but the repulsion and dispersion terms are gathered 
into a Lennard-Jones potential: 
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Finally, the electrostatic interactions are calculated 
with the usual Coulomb potential: 
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Each atom type i is therefore characterized by three 
parameters only, e^, cfi and the parameters of the 
Lennard-Jones potential are then derived using the 
Lorentz-Berthelot mixing rules: 


As in many other fields, most of the classical MD sim¬ 
ulations of RTILs involve non-polarizable force fields and 
they generally include all the atoms from the molecules 
as interaction sites. Since ionic liquids are made of func¬ 
tional groups which are present in many conventional or¬ 
ganic molecules, most of the force fields have been pa¬ 
rameterized following the existing families, such as AM¬ 
BER or OPLS, in which the total interaction potential is 
composed of bonded and non-bonded terms. 


-^total — -^bonds T -^angles T -^dihedrals T -S'non—bonded (1) 

where the bonded terms Ebonds, ^angles and ^dihedrals 


(Jij = ^(cTi + aj) (9) 

The most famous force field for RTILs is the CL&pl^, 
which was built based on the OPLS functional form. 
It was parameterized for a large set of ionic liq¬ 
uids compounds, including the imidazolium, N-alkyl- 
pyridinium, tetraalkylammonium, N,N-pyrrohdinium 
and tetralalkylphosphonium families for cations, and the 
chloride, bromide, triflate, bis(sulfonyl)imide, alkylsul- 
fate, alkylsnlfonate, phosphate, nitrate and dicyanimide 
families for anions.^^^H^ The bonded and Lennard-Jones 
parameters were generally taken form the OPLS sefP^ 
when available and led to accurate results, but some of 







FIG. 2: Snapshots of simulation boxes, using a coloring code 
which identifies the charged (red) and nonpolar (green) do¬ 
mains of RTILs. Left: BMIM-PFe, right: HMIM-PFe. Re¬ 
produced with permission from ref.“ Copyright 2006 Amer¬ 
ican Chemical Society. 


them had to be determined or reparameterized for con¬ 
sistency. This was generally done by reproducing the 
molecular geometries and the torsion energy profiles cal¬ 
culated for isolated molecules using quantum chemistry 
methods (at the MP2 level). The partial charges were 
optimized to reproduce the electrostatic field generated 
by the molecule, using the CHelpG method.!^ 

The first simulation studies of RTILs have focussed on 
the understanding of their struct nr They showed 

that, like in conventional molten salts, it is dominated 
by a combination of short-range repulsion and Coulomb 
ordering effects. Around a given ion, the first neighbour 
shell consists entirely of oppositely charged species; this 
ordering automatically transfers up to several neighbour 
shells. The short-range structure is somewhat more com¬ 
plex: In particular, the three-dimensional arrangement 
of the anions around the cationic molecules strongly de¬ 
pends on the molecular shape of the latter. For example, 
weak hydrogen bonds form between the protons of imi- 
dazoli nm rin gs and the most electronegative atoms of the 
anions 

But the structure of RTILs is also characterized by 
strong intermediate range ordering, which can be de¬ 
tected by the presence of a pre-peak or a shoulder in 
the low wavevector part of X-ray and neutron diffrac- 
tograms. Ribeiro et al showed that both the position 
and the intensity of this prepeak are sensitive to the 
length of alkyl chains of imidazolium cations.^^The pres¬ 
ence of an additional ordering was then clearly character¬ 
ized by Canongia Lopes and Padua. By simulating vari¬ 
ous RTILs using the CL&P and distiguishing the highly 
charged, “polar” regions of the ions from the “unpolar” 
ones, they have put in evidence the presence of nanos- 
tructured domains .1^ Their famous pictures of RTILs are 
reproduced on figure for two different liquids, namely 
the l-butyl-3-methylimidazolium hexafiuorophosphate 
(BMIM-PFe) and the l-hexyl-3-methylimidazolim hex¬ 
afiuorophosphate (HMIM-PFe). This important finding 
is now commonly used to interpret the solvation prop- 
erti es of RTILs towards many spe cies, s uch as nitrate 
ion^HHHor acidic (SO 2 , CO 2 ) gases.It was also ex¬ 


tended to interpret the structure of confined ionic liq¬ 
uids as probed in Atomic Force Microscopy or Surface 
Force Apparatus experiment^^^M^ and in MD simula¬ 
tion^^. The existence of this intermediate range ordering 
has a strong implication on the lubrication properties of 
RTILsP 


III. POLARIZABLE ALL ATOM FORCE FIELDS 



FIG. 3: Gomparison of the simulated and experimental dif¬ 
fusion coefficients for a series of ionic liquids, using po¬ 
larizable (Pol-MD) and non-polarizable (Glass-MD) force 
fields. The corresponding data are respectively extracted from 
referenceJ^ and®. 


Although the CL&P (and related) force fields have led 
to a deep understanding of the structure of RTILs, they 
generally fail to predict the transport properties. Indeed, 
they yield viscosities which are too high by one order 
of magnitude on average, while diffusion coefficients and 
conductivities are underestimated accordingly. This dis¬ 
crepancy is shown on figure where the simulated diffu¬ 
sion coefficients are compared to the experimental ones 
(as measured by pulse field gradient NMR) for a series of 
common ionic liquids (data extracted from referenc^^. 
On this plot, the red line corresponds to an underesti¬ 
mation by a factor of 10. We observe that not only are 
all the values largely underestimated, but also the data 
are largely dispersed and there is not a clear correlation 
between them. This means that even comparing the sim¬ 
ulated diffusivities for two ionic liquids may not give a 
qualitatively meaningful result, and that any dynamic 
information extracted from such a simulation should be 
taken with caution. 

The solution to overcome this difficulty was to include 
polarization effects explicitly in the force field. Indeed, 
since the pioneering study by Yan et al^ there have 
been numerous examples showing an increase of the dif- 
fusivity when using polarizable force fields. Nevertheless, 
these comparisons are often difficult to interpret because 
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FIG. 4: Radial distribution functions between the chloride an¬ 
ions and the protons from the EMIM^ cations in the EMIM- 
AlCU RTIL at 298 K. Top: AIMD results (labelled “CPMD” 
because they were obtained using the Car-Parrinello method); 
bottom: classical MD results with a polarizable force field 
(PIM). Left: Protons from the imidazolium ring. Right: Pro¬ 
tons from the alkyl chains. Reproduced with permission from 
ref.l^. 


one cannot generally add a polarization term on top of 
an existing force field. It is necessary to re-parameterize 
part or all of the other terms. In a recent work, we have 
focussed on the mixtures of l-ethyl-3-methylimidazolim 
chloride (EMIM-Cl) with aluminium chloride (AlCls).^ 
At equimolar proportions, all the chloride ions are linked 
with aluminium ones through very strong ionic bonds, 
thus forming the ionic liquid EMIM-AICI 4 . To simulate 
it, we have used the CL&P force field for all the interac¬ 
tions (both bonded and non-bonded) involving EMIM+ 
cations only. The Polarizable Ion Model (PIM)P^ was 
used for all the other interactions. In this model the re¬ 
pulsion and dispersion terms are described as 
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where the and Cg are the dipole-dipole and dipole- 
quadrupole dispersion coefficients; fn are Tang-Toennies 
dispersion damping functionJ^ describing the short- 
range penetration correction to the asymptotic multipole 
expansion of dispersion. These functions take the form 
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The most important changes concern the electrostatic 
term, which now includes charge-dipole and dipole-dipole 
interaction: 


pPiM 

-^electrostatic 


77 _i_ 77PIM 

-T^^Coulomb I -T^^polarization 


47reo 


EE 

I j>i 


(<m_ 

V nj 


( 12 ) 


+g^\r’^^)qjTYHi - 

+ E5V^? 


where /jLi is the induced dipole of ion i, while Ti and T 2 
are the charge-dipole and dipole-dipole interaction ten¬ 
sors defined by 
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and ai is the polarizability of ion i, which is assumed 
to be isotropic. The last term in equation accounts 
for the energy cost of deforming the electron densities of 
the ions to create the induced dipoles. Note that Tang- 
Toennies functions are also included to account for the 
short-range damping effects on the charge-dipole interac¬ 
tions: 


g^\nj) = 1 - ^ (15) 

k=0 

The parameters in the polarizable potentials are then 
optimized by matching the dipole and forces obtained us¬ 
ing the potential to the ones calculated by DFT on a se¬ 
ries o f repr esentative configurations of liquid chloroalumi- 
nates,^^^^^ following the method described in referenceJ^ 
ancP^. The force field was then validated in two steps: 
Firstly, the structure was compared to that calculated 
with an AIMD simulation. Various radial distribution 
functions provided by the two methods (the AIMD one 
is labelled “CPMD” since it was obtained using the Car- 
Parrinello method) are shown on figure An excellent 
agreement is observed for all the atoms. We note that 
the agreement is much better than the one obtained by 
Youngs et al. who have tried to fit a non-polarizable force 
field for the ionic liquid dimethylimidazolium chloride 
using a similar force-fitting technique.!^ Although they 
were able to reproduce the peak positions better than 
with previous classical MD potentials, the intensities of 
the first peaks where largely overestimated (especially in 
the case of the ring protons). The over-structuring ob¬ 
tained by these authors could partly be cancelled by mul¬ 
tiplying all the €i parameters by a factor of 2 , at the price 
of a worse reproduction of the initial set of DFT forces. 
Taking into account anion polarization effects leads to 
strong improvements in the description of the structure 
of the RTILs. 

















Secondly, the dynamical properties determined with 
the PIM force field for EMIM-AICI 4 were compared to 
the available experimental data. An excellent agreement 
was found for the diffusion coefficients, the viscosity and 
the electrical conductivity,!^ confirming that the polar¬ 
ization term is compulsory for studying transport prop¬ 
erties of RTILs. This is remarkable since the force field 
was parametrized using DFT calculations only, no ex¬ 
perimental information was included in the fitting pro¬ 
cedure. Another interesting aspect of our PIM is that it 
was made consistent with our previous work on molten 
salts,!^ i.e. the chloroaluminate anion is represented by 
A13+ and Cr ions, which are bound together only with 
by strong electrostatic interactions. It could therefore be 
used for example for studying non-stoichiometric mix¬ 
tures of EMIM-Cl and AICI 3 , something which would 
not be possible with a conventional non-polarizable force 
field. 


Polarizable force fields have also been devel- 
opped for a series of RTILs by Borodin, containing 
l-methyl-3-alkyhmidazohum, 1-alky 1-2-methyl-3- 

alkylimidazolium, N-methyl-N-alkylpyrrolidinium, 

N-alkylpyridinium, N-alkyl-N-alkylpiperidinium, 

N-alkyl-N-alkylmorpholinium, tetraalkylammo- 

nium, tetraalkylphosphonium, N-methyl-N- 

oligoetherpyrrolidinium cations and BF 4 , CF 3 BF 3 , 
CH 3 BF 3 , CF 3 SO 3 , PFg, dicyanamide, tri- 
cyanomethanide, tetracyanoborate, bis(trifluoromethane 
sulfonyl)imide (Ntf^ or TFSI“), bis(fiuorosulfonyl)imide 
(FSI“) and nitrate anions.!^ The bonded terms have 
a similar analytical form as in the CL&P while the 
non-bonded terms resemble the PIM ones. For the 
latter, the only difference is that the Tang-Toennies 
functions are replaced by a strongly repulsive term at 
(very) short-range for the dispersion, and through the 
use of Thole screening function^^ for the polariza¬ 
tion. This polarizable force field was parameterized by 
combining quantum chemistry data (for the bonded 
and electrostatic terms) and experimental results for 
the repulsion and dispersion terms. For the latter, 
in addition to the densities, the diffusion coefficients 
were used in the fitting procedure. As shown in figure 
they are very well reproduced by the force field: 
only small discrepancies are observed for the lowest 
values. Such an agreement is of course partly due to the 
parameterization procedure, but it is likely that it would 
not be possible to obtain it using a non-polarizable force 
field. Another example of excellent reproduction of the 
dynamical properties was obtained by Choi et al for the 
BMIM-BF 4 ionic liquid. Their polarizable interaction 
potential was entirely parametrized from ab initio 
calculations, using the symmetry-adapted perturbation 
theory (SAPT). The diffusion coefficients and electrical 
conductivities calculated with this potential are in 
close agreement with the experimental data at several 
temperatures.!^ 


IV. ON THE USE OF REDUCED CHARGES IN 
NON POLARIZABLE FORCE FIELDS 

It is worth noting that in recent years many papers 
have discussed the opportunity to use non-polarizable 
force fields, but with reduced charges. Indeed, although 
the atomic charge is not a well-defined quantity in quan¬ 
tum chemistry or DFT calculations, it is possible to de¬ 
rive some values using a variety of approaches. The most 
popular ones are the CHelpG,!^ RESP,!^ Bader analy- 
si^^ and BlochP^ methods. All these methods, when 
used on ionic liquids with no constraint on the total 
charge of the ions, yield values smaller than the formal 
charge.!^^M^ There are two interpretations: This may be 
either due to polarization or to charge transfer effects. 
Discriminating between these effects is a tricky issue, and 
it is rendered even more difficult by the fact that all the 
charge derivation methods provide different values for the 
same configuration and that even a given method yields 
different sets of charges depending on the conditions of 
the calculation (i.e. choice of the grid points for calculat¬ 
ing the electrostatic potential, different conformations, 
gas-phase or condensed phase calculation, etc). Never¬ 
theless the success of the PIM for a large variety of inor¬ 
ganic ionic material^ and of polarizable force fields for 
RTILs seem to indicate that charge transfer is negligible 
in these systems. An advantage of the SAPT approach 
by Choi et al is that the induction energy, which con¬ 
tains both the polarization and the charge transfer terms, 
is calculated with ab initio techniques. In their study, 
they concluded that polarization is the dominant energy 
based on the fact that it is extremely well-reproduced 
by their force field which includes polarization only (the 
ab initio calculation was performed on 1300 BMIM-BF 4 
complexes extracted from the liquid).!^ 

Using reduced charges should therefore be seen as an 
effective way to account for polarization effects only. A 
rationale for this approach was provided by Leontyev 
and Stuchebrukhov, who explain the need for charge¬ 
scaling in non-polarizable force fields by a neglect of the 
electronic solvation energy.!^ They therefore suggest to 
screen the Coulomb interactions by a factor 1 — l/eez, 
where Cqi is the electronic (high-frequency) dielectric con¬ 
stant, which can be determined from experimental mea¬ 
sures of the refractive inde:x!^ or from a direct compu¬ 
tation of the molecular polarizabilities using DFT.!® 
order to check how reliable this simplified approach is, 
Schroder performed an extensive comparison of the two 
methods. In his study, he progressively “switched on” the 
polarization effect (using Drude oscillators instead of the 
explicit induced dipoles discussed here, but this should 
not impact much the conclusion J®) , either by increasing 
the polarizability or by scaling down the charges with a 
scaling factor / ranging from 0 to unity (i.e. the sim¬ 
ulation with / = 0 corresponds to the non-polarizable 
force field, while / = 1 corresponds to the normal polar¬ 
izable force field or to the non-polarizable force field with 
scaled charges).!® He showed that although the reduced 


charge force field was not able to reproduce the correct 
dipole distribution, leading to substantial deviations for 
the mean rotational relaxation time, the diffusion coeffi¬ 
cients and the electrical conductivity were qualitatively 
correct. In conclusion, the reduced charge method can 
be considered appropriate if one is interested in recov¬ 
ering the correct order of magnitude for the dynamical 
properties, but it is important to keep in mind that the 
local relaxation mechanisms may not be accurate. 

As discussed in the introduction, an important aspect 
of a simulation is its computational cost. Including ex¬ 
plicit dipoles increases the simulation times by a fac¬ 
tor of 5 to 10 for a given number of atoms, which ex¬ 
plains why they are often discarded. Nevertheless, we 
observe on figure that this order of magnitude in the 
cost does not seem to be reflected in practice. Of course, 
the data shown here is far from being exhaustive and 
there are many possible biases (for example, it is possible 
that the groups who have access to larger computational 
ressources will tend to use a more expensive method for 
tackling a given problem), but the expected difference of 
one order of magnitude is not clearly observed. Indeed, 
although larger systems are generally simulated by non- 
polarizable force fields, the longer trajectories have been 
accumulated using polarizable ones. These observations 
may be explained by the fact that the polarizable force 
fields are generally used for calculating transport prop¬ 
erties while the non-polarizable ones focus on the struc¬ 
ture. Note also that in classical MD the simulation time 
scales as where N is the number of interacting sites, 
while extending a trajectory over time leads to a linear 
increase of the simulation time. As a consequence, it is 
more efficient to do a longer simulation time on a smaller 
simulation cell for studying the transport properties with 
polarizable force fields. Nevertheless figure shows that 
the continuous increase in the computational ressources 
available, together with the efficient parallelization of MD 
codes, has allowed polarizable force fields to become a 
viable alternative for studying the physical properties of 
RTILs. 


V. COARSE-GRAINED FORCE FIELDS 

Coarse-grained force fields are used when the compu¬ 
tational cost needs to be strongly reduced. For example, 
many problems in material science in which ionic liquids 
are used involve complex interfaces. In energy storage 
devices such as supercapacitors or batteries, RTILs play 
the role of the electrolyte, in which the transport of the 
charge occurs via the migra tion of ionic species from one 
electrode to the other.I^SSl Understanding how they op¬ 
erate then necessitates very large simulation cells with 
explicit electrodes. This limits the number of systems 
that can be studied and prevents systematic comparisons. 
Here we will focus on the topic of super capacitors, which 
has attracted a lot of attention from the RTIL modelling 
community in the last years. Indeed, the energy stor¬ 


age does not imply any (electro)chemical reaction, but 
rather the reversible adsorption of the ions on porous 
carbon electrodes, so that this problem can be tackled 
using classical molecular simulation methods. A first un¬ 
derstanding of the adsorption behavior of RTILs on elec¬ 
trodes has been provided by analytical theories .1^ Then 
the first simulations have been performed using simpli¬ 
fied models of RTILs. For example, Fedorov et al. de¬ 
vised some generic rules on the behavior of RTILs at the 
surface of planar electrodes, with an emphasis on the ef¬ 
fects of the ionic sizes, of the presence/absence of neutral 
groups (represen ting the alkyl chains) on the cations and 
of the voltage.l^^^^ Nevertheless it is worth underlining 
that these simplified models cannot really be considered 
as coarse-grained force fields because there is no map¬ 
ping between them and real ionic liquids, so that it is 
not possible to use them directly for the study of a spe¬ 
cific system. 



FIG. 5: Coarse-graining of the EMIM^-BF 4 ionic liquid, a) 
Electron density isosurface mapped with the corresponding 
electrostatic potential, b) All-atom representation, c) Coarse¬ 
grained representation. 

The two main challenges for simulating realistic su¬ 
percapacitors are i) the size of the system and ii) the 
handling of constant potential cond itions . Indeed, they 
involve nanoporous carbon materials) ^^ * ^^ * which can only 
be modelled with complex structures. There are several 
methods for maintaining those at constant potential.!^ 
In our model,I^^M^ the local charges on the atoms of the 
electrode vary dynamically in response to the electrical 


potential caused by the ions and molecules in the elec¬ 
trolyte. They are therefore calculated at each time step 
of the simulation using a self-consistent approach simi¬ 
lar to the one used for calculating the induced dipoles in 
polarizable force fields, so that the computational cost 
associated with such simulations is very high. For these 
reasons, simulations at constant applied potential and us¬ 
ing all atom force fields have mostly been limited to elec¬ 
trodes with simpler geometries and the computational 
cost is often reduced by using multiple time step algo¬ 
rithms, i.e. the values of the electrode ch arge s are not 
calculated at each step of the simulation.^^^^i^ 

In our recent work we have focussed on supercapacitors 
involving carbide-derived carbon electrodes and RTIL 
electrolytes, which present very promising experimental 
performance s 1^^1 1511 The carbon structures, obtained by 
Palmer et alp^sne made of around 4000 atoms and have 
pore size distri butio ns and accessible surfaces similar to 
the real devicesDue to the large number of electrode 
atoms, it was not possible to use an all atom force field 
for the RTIL, so that we resorted to a coarse-grained 
reprensentation of the ions instead. Experimental works 
have shown that the most important effect on the capac¬ 
itance is d ue to th e matching between the pore size and 
the ion siz d^^^ * ^^ ^, so this approximation should not af¬ 
fect the results from the qualitative point of view. The 
possibilities to describe a RTIL using a coarse-grained 
approach are very large. Prior to the choice of the func¬ 
tional form for the non-bonded terms (and to the deter¬ 
mination of the corresponding parameters), it is neces¬ 
sary to establish which atoms will be put together in a 
grain and whether the “bonds” should be held rigid or 
treated with bonded terms similar to those of equations 
[^[^and|^ In a first step, simple RTILs based on BMIM+ 
and EMIM+ cations and on PEg and BE 4 anions were 
considered. As can be seen on the panel a) of figure 
the electrostatic potential (extracted from a DET calcu¬ 
lation) of these molecules is rather smooth, so that it is 
possible to pass from the all-atom representation of panel 
b) to the coarse-grained one of panel c). Such a coarse- 
grai ning approach was proposed for BMIM-PEe by Roy 
et In their representation, the BMIM+ cation is 

represented by three interaction sites (one for the imida- 
zolium ring, one for the methyl group and one for the 
butyl group) while the PE^ anion is replaced by a sin¬ 
gle site. The cation is kept rigid during the simulation, 
its geometry is therefore chosen according to its aver¬ 
age conformation in the liquid. The non-bonded interac¬ 
tion are calculated using Lennard-Jones (equation]^ and 
Coulomb (equation ^ potentials. The corresponding pa¬ 
rameters (e^, ai and g^) do not have an atomic character 
anymore, so that they are less transferable. Starting from 
an initial guess based on an all-atom force field, they 
are chosen in order to match as closely as possible the 
experimental density and diffusion coefficients, so that 
they are purely empirical. In a first step, unit charges 
were attributed to the ions.^^^ As for all atom non po¬ 
larizable force field, this led to understimated transport 


properties. The use of reduced charge^^^ then corrected 
this discrepancy. In the case of a coarse-grained model, 
this choice is fully justified since there should be a strong 
screening from the whole grain associated to each charge. 
This coarse-grained force field was t hen validated for in¬ 
terfacial properties (s urfac e tension]P^ and extended to 
EMIM+ and BE 4 ionP^l. 


FIG. 6: Instantaneous structure of the BMIM-PFe ionic liquid 
adsorbed in a nanoporous carbon electrode held at 0.5 V. Left: 
Position of the atoms (turquoise: carbon, green: PF^, red: 
BMIM^). Right: Local charge of the various carbon atoms 
(green: negative charge, red: positive charge and yellow: null 
charge); the atoms of the ionic liquids are shown in gray. 

Using this potential, we have been able to obtain much 
larger capacitances in nanoporous carbons than with pla¬ 
nar elect rodes, in agreement with the experimental re- 
sults.ll^ We have shown that the electrode is wetted by 
the electrolyte at null potential and that the charging 
process involves the exchange of ions with the bulk elec¬ 
trolyte without substantially changing the volume of liq¬ 
uid inside the electrode. This exchange is accompanied 
by a partial decrease of the coordination number of the 
ions rendered possible by the charge compensation by the 
electrode, which underlines the importance of the elec¬ 
trostatic interactions on the charge storage mechanism. 
The screening of the Coulomb interaction by the metal 
at non-zero voltages is even at the origin of the forma¬ 
tion of a “superionic” state,^^^^^^^ in which two cations 
or two anions are contiguous as can be seen on figure a 
situation which is extremely unfavorable in a bulk RTIL. 
Then we demonstrated the key role of the local str ucture. 
In qualitative agreement with recent NMR studies,^!!^®^ 
our simulations could put in evidence that the adsorption 
of the ions occurs on vari ous sites in a way which depends 
on the applied potential.^^^ Four different types of sites 
were proposed in the case of CDCs: Depending on the 
number of carbon atoms in the vicinity of the ions, we 
have distinguished edge sites (concave curvature), plane 
sites (graphene sheet-like structure), hollow sites (convex 
curvature) and pocket sites (inside a subnanometre car¬ 
bon pore). More recently, we could exploit the fact that 
these coarse-grained force fields were adjusted in order to 
reproduce well the diffusion coefficients for stud ying the 
charging dynamics of these super capacitors Charac¬ 
teristic charging times of approximately 1 second were 











determined by extrapolating the MD results to a macro¬ 
scopic electrode, here again in good agreement with the 
available experimental data. 

Although very efficient, these coarse-grained force 
fields suffer from a lack of transferability. It is probably 
necessary to reparameterize them and to redefine the ge¬ 
ometry of the ions when passing from a RTIL to another. 
If more complex ionic liquids are simulated in the future, 
it will probably be necessary to use the effectiv e force 
coarse-graining approach proposed by Voth et 
It is rather similar to the force-fitting approach that we 
have used for polarizable force field. The idea is to de¬ 
termine effective pairwise forces between coarse-grained 
sites by averaging over the atomistic forces between the 
corresponding atomic groups in configurations extracted 
from all-atom MD simulations (in principle it would also 
be possible to use AIMD). The size of the grains is smaller 
than in the models proposed by Roy et a/., and the poten¬ 
tial includes bonded interaction terms, so that it should 
be possible to keep similar parameters for a wide family 
of ionic liquids. This approach is therefore very promis¬ 
ing as a next step for the study of super capacitors. 

VI. CONCLUSIONS 

Since its burgeoning, the field of ionic liquids has taken 
much profit from the use of molecular dynamics. On 
the one hand, AIMD provides an unique framework for 
studying the vibrational properties and the mechanisms 
of chemical reactions. On the other hand, classical MD 
allows to reach the time and length scales necessary for a 
deep understanding of their structural, thermodynamic 
and transport properties, both in the bulk and at inter¬ 
faces. In the past ten years, such simulations have pro¬ 
vided a quantitative understanding of the structure of the 
ionic liquids, showing the formation of nanosegregated 
polar/apolar domains, which has strong consequences on 
the solvation properties of these media. The determi¬ 
nation of transport properties has proven more difficult: 
conventional, non-polarizable force fields largely underes¬ 
timate the diffusion coefficients and polarization effects 
need to be introduced to recover the correct dynamical 
behavior. In recent years, two approaches have been pro¬ 
posed. The simplest one is to use rescaled charges for the 
ions, while the other one consists in introducing explicit 


dipoles (either using Drude oscillators or point dipoles). 
Although both approaches provide enhanced dynamics, 
the latter is more accurate. In addition, the correspond¬ 
ing parameters are more transferable from one compound 
to another, so that this approach should be preferred for 
systematic studies. Indeed, an ambitious but feasible ob¬ 
jective for future studies should be the prediction of the 
physico-chemical properties of RTILs prior to their syn¬ 
thesis, in order to target the most adapted one for a given 
application. 

Applications, however, often use ionic liquids under 
special conditions. Indeed, they are widely used as elec¬ 
trolytes in energy storage devices, and it is then necessary 
to understand their behavior at complex interfaces. The 
computational cost associated to heterogeneous systems 
is usually very large, which implies that these studies are 
currently tackled using coarse-grained force fields. The 
loss of atomic details implies that comparisons between 
specific RTILs become more difficult, but these simula¬ 
tions allow for the understanding of problems which are 
otherwise difficult to address by experiments only. For 
example, the field of supercapacitors has recently ben- 
efitted greatly from molecular simulations, which have 
provided interpretations for the changes in the electric¬ 
ity storage ability depending on the geometry of the elec¬ 
trodes based on the structure of the ionic liquids. The 
dynamical aspects are more difficult to handle since even 
with coarse-grained force fields the computational costs 
remain high. There are many other fields in materials 
science in which molecular simulations of RTILs at inter¬ 
faces should be able to provide ve ry u seful input for op- 
timizing device s: e lectroactuation,^^^ lubrication, b^QF2i ] 
Li-ion batteries 
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